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2.  Objectives 

This  project  is  concerned  with  the  development  of  finite  element  procedures  for 
effectively  performing  metal  forming  simulations  in  an  automated  analysis  envi¬ 
ronment.  Due  to  their  ability  to  efficiently  provide  results  to  the  required  levels  of 
accuracy,  the  desired  formulation  should  support  hp-adaptive  analysis  techniques. 
These  requirements  present  a  technical  challenge  due  to  the  inability  of  standard 
finite  element  formulations  to  satisfy  the  Babuska-Brezzi  conditions  for  general 
combinations  of  displacement  (velocity)  and  pressure  interpolations  needed  for  hp- 
adaptivity  in  the  presence  of  incompressibility  constraints.  To  address  this  problem 
we  have  investigated  stabilized  finite  element  formulations  for  this  class  of  prob¬ 
lem.  The  second  area  considered  was  the  development  of  a  geometry-based  simu¬ 
lation  framework  capable  of  supporting  automated  hp-adaptive  technologies  using 
such  element  formulations  on  problems. 

3.  Stabilized  Finite  Element  Formulation 

Typical  metal  forming  processes  involve  large,  nearly  incompressible,  deforma¬ 
tions.  In  order  to  handle  the  incompressibility,  mixed  formulations,  where  the  dis¬ 
placements  and  pressures  are  interpolated  separately,  have  proved  to  be  effective. 
However,  it  has  been  shown  that  these  mixed  formulations  are  not,  in  general,  sta¬ 
ble  unless  they  satisfy  a  certain  stability  criterion,  the  so-called  Babuska-Brezzi 
condition.  This  criterion  puts  restrictions  on  the  relationship  between  the  displace- 
•  xnent.  and  pressure  finite  element  interpolation  functions.  This  restriction  limits 
computational  efficiency  since  the  order  of  interpolation  is  defined  by  this  stability 
condition  rather  than  by  numerical  accuracy,  and  thus  makes  p-adaptivity  all  but 
impossible.  Because  metal  forming  processes  are  generally  three-dimensional  in 
nature  and  involve  very  large  deformations,  adaptivity  is  critical  to  obtaining  accu¬ 
rate  results  efficiently.  A  combination  of  h-  and  p-  adaptivity  is  most  desirable  for 
gaining  optimal  efficiency.  In  this  work,  a  stabilized  finite  element  formulation  for 
large  deformation  processes  which  eliminates  the  need  to  satisfy  the  Babuska- 
Brezzi  condition  is  presented.  The  formulation  is  based  on  the  stabilized  formula¬ 
tions  developed  for  linear  problems  by  Hughes  et  al.  Examples  involving  large 
deformation,  three-dimensional,  hyperelasticity  with  linear  interpolations  for  both 
the  displacement  and  pressure  are  presented  to  demonstrate  the  algorithm.  The 
results  show  that  the  method  can  work  for  large  deformation,  nonlinear,  three- 
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dimensional  problems,  and  is  therefore  a  promising  technique  for  effectively  mod¬ 
eling  metal  forming  processes. 

Given  a  continuum  body  iP  we  introduce  the  deformation  cp.P  — >  A  and  use  fee 
notation  B  =  <p  (?)  c  9T  for  fee  reference  configuration  and  B  =  <p (P)  c  A 
for  the  current  configuration.  A  point  in  the  reference  configuration  X  is  mapped 
into  a  point  *  =  <p(X)  =  X  +  u(X)  in  the  current  configuration  where  u  denotes 
the  displacement  field.  We  define  the  deformation  gradient 

p  =  Vcp(.X’)  =  1  +  Vn  and  use  the  notation  C  -  F  F  for  the  right  Cauchy- 
Green  tensor.  The  symbols  V[«]  and  Div[«]  denote  the  gradient  and  divergence, 
respectively  of  a  quantity  [•]  with  respect  to  the  coordinates  X  e  B0.  Assuming 
the  volumetric  deformation  is  purely  elastic,  we  can  write  p  =  kU'(J)  ,  where  p  is 
the  mean  stress,  k  is  the  bulk  modulus,  and  U  is  the  volumetnc  part  of  the  stored 
energy  depending  only  on  fee  Jacobian  J  of  the  deformation  Gradient  F .  The  sec¬ 
ond  Piola-Kirchhoff  stress  S  can  be  decomposed  into  S  =  JpC  +S  where  the 
first  term  is  the  volumetric  part  and  S  is  the  remainder.  In  fee  following  we  present 
a  stabilized  mixed  displacement-pressure  formulation  in  matenal  quantities 
defined  in  the  reference  configuration.  The  boundary  value  problem  for  equilib¬ 
rium  in  the  absence  of  body  forces  is  given  as:  Find  a  displacement  field  u  such 

that  ■ 

-Div  [FS]  =  0  in  B0 

[FS]n0  =  g0  on  ToN  .  (D 

u  =  u0  onToD 

n  is  the  exterior  unit  normal  on  Ta ,  the  boundary  of  B0 ,  and  g0  is  a  prescribed 
traction  load  on  part  of  fee  boundary  VoN .  We  consider  prescribed  dmchlet  bound¬ 
ary  conditions  u,  on  ToD.  ToD  and  ToN  describe  the  complete  boundaij  of  B 
The  stabilized  mixed  finite  element  formulation  can  now  be  derived  from  the 
strong  form  Eqn.m  by  multiplying  wife  the  perturbed  weighting  function 
u  +  (a/i2)/(2p)F  Jr  v| ,  and  integrating  over  fee  domain 

\  Div[F5]  •  udV  +  X  ^  j  Div£FSl '  F_Tv?dy  =  0  (2) 

s.  e  =  1  *.  _T 

where  fee  perturbation  term  (ah2)/ (2p)F  Vp  is  applied  element  wise,  a  is  a 
stabilization  parameter,  h  is  a  mesh  size  parameter,  and  \x  is  shear  stiffness.  The 
first  term  in  Eqn.  (2)  is  integrated  by  parts  as  usual.  We  will  focus  our  attention  on 
the  derivation  of  the  stabilization  term.  Using  fee  additive  decomposition  of  the 
second  Piola-Kirchhoff  stress  we  obtain 

Div[FS]  =  Divtp/FC-1]  +Div[F5]  (3) 

=  7F-TVp  +Div[FS] 

The  second  step  in  the  formulation  was  achieved  by  applying  the  Piola  Identity 
Di w[JF~T]  =  0 .  Introducing  Eqn.  (3)  into  Eqn.  (2)  and  introducing  the  pressure 


p  =  K.lf(J(u))  as  an  independent  variable,  the  stabilized  mixed_weak  formula¬ 
tion  of  Eqn.  (1)  is  then:  Find  ( u,p)eVxQ  such  that  for  all  {u,p)eVxQ 


|5:[FrVw]dV+  J  JpC  l:[FTVu]dV  =  Lext{u) 
B„  Bo 


In  the  case  when  only  linear  shape  functions  are  used,  the  divergence  term  in  Eqn 
(4)  vanishes.  We  consider  a  hyperelastic  material  model  with  an  additive  stored 
energy  function  W  =  k  U(J)  +  W ( C )  so  that 


2k  U\J)  • 


dJ  ,  JW{C) 
dC  dC 


We  choose  a  Neo-Hookean  material  model  where 
W{C)  =  1/2  -  M-[J_2/3trC  -  3] 


=  JpC  l+s, 

U(J )  =  1/2 -(/-l)2 


(5) 

and 


More  complete  details  on  the  formulation  and  numerical  results  showing  the  supe¬ 
riority  of  the  formulation  are  presented  in  references  [18,19]. 


4.  Geometry-Based  Object-Oriented  Simulation  Framework 

To  date  consideration  of  oriented  programming  in  simulation  software  has  focused 
on  flexible  structures  with  code  reuse,  application  of  symbolic  computing  operat¬ 
ing  in  parallel,  linking  with  design  processes  and  supporting  interacting  multiphys¬ 
ics0  simulations.  Building  on  these  efforts  and  the  needs  of  adaptive  simulation 
technologies  we  have  constructed  a  geometry-based  simulation  frameworks  that 
supports  parallel  adaptive  simulation  capabilities.  This  system,  referred  to  as  Trel¬ 
lis  is  based  on  [2,4]:  , 

•  A  set  of  geometry-based  structures  which  can  support;  (i)  the  direct  linkage  ^ 
with  company  CAD  information,  (ii)  all  forms  of  adaptivity  without  introducing 
geometric  approximation  errors  [8],  and  (iii)  the  high  level  integration  of  multi¬ 
scale,  multi-physics  analysis  methodologies. 

•  A  careful  decomposition  of  the  geometry,  physics,  mathematical  model,  discret¬ 
ization  and  numerical  methods  into  interacting  classes.  These  structures  support 
a  variety  of  equation  discretization  methods.  Both  finite  element  [18,19]  and 
partition  of  unity  methods  have  been  implemented  [20]. 

•  Adaptive  control  of  each  step  of  the  simulation  process  from  the  selection  of  the 
mathematical  model,  through  the  model  and  domain  discretization,  to  the  selec¬ 
tion  of  application  of  the  numerical  methods  to  solving  the  discrete  system. 

Conceptually  Trellis  is  built  on  the  view  of  an  analysis  as  a  transformation  between 
three  levels  of  description.  The  highest  level  description  is  that  of  the  physical 
problem  which  is  posed  in  terms  of  physical  objects  interacting  with  their  environ- 


ment.  Since  the  goal  of  the  analysis  is  to  obtain  reliable  estimates  of  the  response 
of  the  system  the  second  level  is  a  mathematical  problem  description  that  intro¬ 
duces  some  level  of  idealization,  which  also  needs  to  be  controlled  to  yield  the 
desired  accuracy.  The  third  level  is  the  numerical  discretization  constructed  from  a 
mathematical  problem  that  involves  another  set  of  idealizations  which  also  need  to 
be  controlled. 

The  structures  used  to  support  the  problem  definition,  the  discretizations  of  the 
model  and  their  interactions  are  central  to  Trellis.  The  two  structures  of  the  geo¬ 
metric  model  and  attributes  are  used  to  house  the  problem  definition.  The  analysis 
discretizations  are  housed  in  the  mesh  structure.  The  final  structure  is  the  field 
structure  which  houses  numerical  solution  results. 

The  geometric  model  representation  is  a  non-manifold  boundary  representation 
based.  The  representation  used  for  a  mesh  is  similar  to  that  used  for  a  geometric 
model  [3]:  a  hierarchy  of  regions,  faces,  edges  and  vertices.  In  addition,  each  mesh 
entity  maintains  a  relation,  called  classification,  to  the  model  entity  that  it  was  cre¬ 
ated  to  partially  represent.  Understanding  how  the  mesh  relates  to  the  geometric 
model  is  critical  for  both  mesh  adaptivity  and  understanding  how  the  solution 
relates  back  to  the  original  problem  description.  The  topological  representation  can 
also  be  used  to  great  advantage  in  performing  adaptive  p-version  analyses  as  poly¬ 
nomial  orders  can  be  directly  assigned  to  the  various  entities  [8,22]. 

A  problem  with  many  “classic”  numerical  analysis  codes  is  that  the  solution  of  an 
analysis  is  given  in  terms  of  the  values  at  a  set  of  discrete  points.  Trellis  eliminates 
this  problem  by  introducing  a  construct  known  as  a  field  which  describes  the  varia¬ 
tion  of  a  tensor  over  one  or  more  entities  in  a  geometric  model.  The  spatial  varia¬ 
tion  of  the  field  is  defined  in  terms  of  interpolations  defined  over  a  discrete 
representation  of  the  geometric  model  entities,  which  can  be  a  mesh. 

The  Trellis  analysis  process  is  a  series  of  transformations  of  the  problem  from  the 
original  mathematical  problem  description  through  to  sets  of  algebraic  equations 
approximately  representing  the  problem.  The  mathematical  problem  description 
level  is  described  by  a  ContinuousSystem  class,  which  contains  the  geometric 
model  and  the  attributes  which  apply  to  that  model,  specified  by  a  particular  case 
node  in  the  attribute  graph.  An  instance  of  a  ContinuousSystem  is  then  transformed 
to  an  instance  of  the  class  DiscreteS y stem  which  represents  the  discretized  version 
of  the  model  and  attributes  and  the  weak  form  of  the  partial  differential  equation 
(PDE).  The  particular  analysis  class  that,  is  used  depends  on  the  selected  weak 
form  of  the  PDE  to  be  solved. 

The  DiscreteSystem  class  represents  the  problem  in  terms  of  contributions  from  a 
set  of  objects  that  live  on  the  discrete  representation  of  the  model.  These  objects 
are  called  SystemContributors.  There  are  three  types  of  SystemContributors.  Stiff- 
nessContributors  contribute  coupling  terms  between  degrees  of  freedom  of  the  sys¬ 
tem,  ForceContributors  contribute  terms  to  the  right  hand  side  vector,  and 
Constraints  set  specific  values  or  constraints  to  given  degrees  of  freedom.  These 
objects  are  created  by  the  Analysis  object  and  correspond  to  an  interpretation  of 
attributes  consistent  with  the  weak  form  that  the  Analysis  implements. 


The  Analysis  class  creates  all  of  the  SystemContributors  and  adds  them  to  an 
instance  of  a  DiscreteS ystem.  The  DiscreteSystem  is  transformed  into  an  Algebra- 
icSystem,  an  Assembler  object.  Multiple  linear  solvers  can  be  used  to  solve  the 
Algebraic  System.  The  most  extensive  capability  included  is  the  Portable  Extensi¬ 
ble  Toolkit  for  Scientific  Computation  (PETSc)  from  Argonne  National  Labora¬ 
tory.  These  procedures  have  the  dual  advantage  of  working  effectively  in  an  object- 
oriented  analysis  framework  and  providing  an  efficient  set  of  linear  algebra  rou¬ 
tines. 
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